library(readr)    
library(readxl)     
library(stringr)    
library(dplyr)      
library(xtable)     
library(writexl)    
library(Matrix)     
library(ivreg)      
library(reshape2)   
library(SQUAREM)

load("~/Salidas/b. MDATA4.RData")
DATA <- DATA %>% filter(EMPRESA != "VIRGIN")

ns <- 10000
set.seed(123)
v <- matrix(rnorm(ns), nrow(DATA), ns, byrow = T)
save(v, file = "c. v_02.RData")

DATA <- DATA %>%
  arrange(t, EMPRESA, serv, posp) %>%
  group_by(t) %>%
  mutate(
    So = 1 - sum(s, na.rm = TRUE)
  ) %>%
  ungroup()

p      <- DATA$p     			
S      <- DATA$s     			
So     <- DATA$So    			
delta0 <- log(S) - log(So)      	
t      <- DATA$t     			
J      <- nrow(DATA) 			
NS     <- ncol(v)    			
mkti   <- match(DATA$t, unique(DATA$t)) 

X <- DATA %>%
  transmute(
    posp = posp
  ) %>% as.matrix()

Xd <- DATA %>%
  transmute(
    posp = posp
  ) %>% as.matrix()

dx <- DATA %>%
  transmute(
    p, 
    XB = paste0('_', j),       
    XF = paste0('_', EMPRESA), 
    XT = paste0('_', ANNO),    
    Xt = paste0('_', t)        
  )

XB <- model.matrix(~ 0 + XB, data = dx)[,-1]
XF <- model.matrix(~ 0 + XF, data = dx)[,-1]
XT <- model.matrix(~ 0 + XT, data = dx)[,-1]
Xt <- model.matrix(~ 0 + Xt, data = dx)[,-1]

Tn <- length(unique(DATA$t))
H  <- matrix(0, J, J)
m  <- unique(DATA$t)

for(tt in 1:Tn) {
  idx <- which(mkti == tt)
  x_f <- XF[idx, , drop = FALSE]
  x_f <- cbind(base_firm = 1 - rowSums(x_f), x_f)
  H[idx, idx] <- x_f %*% t(x_f)
}

id <- c(2, 3, 5)
Z <- DATA %>%
  select(all_of(paste0('z', id))) %>%
  as.matrix()
Z <- cbind(Z, constant = 1, X, XF, Xt)
A <- solve(t(Z) %*% Z, tol = 1e-30)

XPF <- cbind(constant = 1, X, XF, XT)
WPF <- cbind(constant = 1, X, XF, XT)

ZPZ  <- Z %*% A %*% t(Z)
ZPZ2 <- bdiag(ZPZ, ZPZ) %>% as.matrix() 
XZPZ <- solve(t(XPF) %*% ZPZ %*% XPF) %*% t(XPF) %*% ZPZ
WZPZ <- solve(t(WPF) %*% ZPZ %*% WPF) %*% t(WPF) %*% ZPZ

zeta0 <- c(0.1, 0.3)
lbs <- c(-3, 0.001)
ubs <- c(3, 3)
ans1 <- nlminb(start = zeta0, objective = gmm, lower = lbs, upper = ubs, control = list(iter.max = 10000, eval.max = 500, trace = 1))

zeta <- ans1$par
out  <- struct(zeta) 	     
b    <- out$b

zeta1 <- ans1$par
xi    <- out$xi  	     
wo    <- out$wo  	     
om    <- c(out$xi, out$wo)   

A2    <- 0
for(i in 1:nrow(Z)) A2 <- A2 + om[i]^2 * Z[i, ]%*%t(Z[i, ])/(nrow(Z)*3)
for(i in 1:nrow(Z)) A2 <- A2 + om[i + nrow(Z)]^2 * Z[i, ]%*%t(Z[i, ])/(nrow(Z)*3)
A2    <- solve(A2, tol = 1e-30) 

ZPZ   <- Z %*% A2 %*% t(Z)
ZPZ2  <- bdiag(ZPZ, ZPZ) %>% as.matrix()

lbs <- c(-3, 0.001)
ubs <- c(3, 3)
ans <- nlminb(start = zeta1, objective = gmm, lower = lbs, upper = ubs, control = list(iter.max = 10000, eval.max = 500, trace = 1))

zeta  <- ans$par
out   <- struct(zeta)

delta <- out$delta
cm    <- out$cm
si    <- out$si
b     <- out$b; b 
ai    <- -(zeta[2] + zeta[1]*v)
g     <- out$g; g 
xi    <- out$xi   
wo    <- out$wo   

save(list = c('cm', 'delta', 'zeta'), file = "c. Estimout_01_06.RData")

om       <- c(xi, wo)
Z2       <- rbind(Z, Z)
theta    <- c(zeta, b)
eps      <- 1e-13
Kg       <- ncol(Z)
Kb       <- length(theta)
dtheta   <- numeric(Kb)
G   <- PSI <- 0
OM1 <- OM2 <- matrix(0, length(om), Kb)

for(k in 1:Kb) {
  dtheta[k] <- eps
  if(k <= length(zeta)) {
    om1 <- omega2(zeta = zeta - dtheta[1:length(zeta)])[1:nrow(Z2)]
    om2 <- omega2(zeta = zeta + dtheta[1:length(zeta)])[1:nrow(Z2)]
  } else {
    om1 <- omega1(theta[-(1:length(zeta))] - dtheta[(length(zeta) + 1):Kb])[1:nrow(Z2)]
    om2 <- omega1(theta[-(1:length(zeta))] + dtheta[(length(zeta) + 1):Kb])[1:nrow(Z2)]
  }
  dtheta  <- numeric(Kb)
  OM1[, k] <- om1
  OM2[, k] <- om2
}

for(i in 1:nrow(Z2)) {
  print(i) 
  Gi <- matrix(0, Kg, Kb)
  for(k in 1:Kb) {
    g0 <- Z2[i, ] * om[i]
    g1 <- Z2[i, ] * OM1[i, k]
    g2 <- Z2[i, ] * OM2[i, k]
    Gi[, k] <- 0.5*(g1 - g0)/(-eps) + (g2 - g0)*0.5/eps 
  }
  G    <- G + Gi/nrow(Z2) 
  PSI <- PSI + om[i]^2 * Z2[i, ]%*%t(Z2[i, ])/nrow(Z2)
}

W <- A2

A   <- t(G) %*% W %*% G
A   <- solve(A, tol = 1e-20)
B   <- t(G) %*% W %*% PSI %*% W %*% G
Var <- (A %*% B %*% A) / (J - Kb)

se <- sqrt((diag(Var)))
bb <- c(zeta, b)
Kb <- length(bb)
sb <- tail(se, length(bb))
pv <- (1 - pt(abs(bb/sb), df = J - Kb))/2
x  <- cbind(bb, sb, pv)

xr <- x[c(2,1), ]
rownames(xr) <- c('alpha', 'sigma')
print(xtable(xr, digits = 3))

xr <- x[c(3:7), ]
rownames(xr) <- paste0('beta-', c('constante', 'pospago', 'MOVISTAR', 'TIGO', 'PTC'))
print(xtable(xr, digits = 3))

xr <- x[c(8:10), ]
rownames(xr) <- paste0('beta-', c('2023', '2024', '2025'))
print(xtable(xr, digits = 3))



